﻿using System.Numerics;
using WaveSim.Core.Equations;

namespace WaveSim.Core.Methods.ODE
{
    [MethodODE("Adams 4", Description = "The Implicit Adams method of 4th degree")]
    public class ImpAdams : IMethodODE
    {
        private Complex[,] F;
        Complex[] yt, k1, k2, k3, k4;

        private void init(Complex[] y)
        {
            if (F == null)
            {
                F = new Complex[y.Length, 4];
            }
        }

        private void init(int n)
        {
            if (yt == null)
            {
                yt = new Complex[n];
                k1 = new Complex[n];
                k2 = new Complex[n];
                k3 = new Complex[n];
                k4 = new Complex[n];
            }
        }

        public void Solve(double x, double h, Complex[] y, IEquation f, int iteration = 0)
        {
            if (iteration < 3)
            {
                init(y);
                // calculate value using Runge-Kutta
                init(y.Length);
                for (int i = 0; i < y.Length; i++)
                {
                    k1[i] = h * f.Get(i, x, y);
                    yt[i] = y[i] + 0.5 * k1[i];
                }
                for (int i = 0; i < y.Length; i++)
                {
                    k2[i] = h * f.Get(i, x + h * 0.5, yt);
                    yt[i] = y[i] + 0.5 * k2[i];
                }
                for (int i = 0; i < y.Length; i++)
                {
                    k3[i] = h * f.Get(i, x + h * 0.5, yt);
                    yt[i] = y[i] + k3[i];
                }
                for (int i = 0; i < y.Length; i++)
                {
                    k4[i] = h * f.Get(i, x + h, yt);
                }
                for (int i = 0; i < y.Length; i++)
                {
                    y[i] += (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6;
                }

                for (int i = 0; i < y.Length; i++)
                {
                    F[i, iteration] = f.Get(i, x + h, y);
                }
            }
            else
            {
                var copy = (Complex[])y.Clone();

                for (int i = 0; i < y.Length; i++)
                {
                    y[i] += (55 * F[i, 3] - 59 * F[i, 2] + 37 * F[i, 1] - 9 * F[i, 0]) * h / 24;
                }

                for (int i = 0; i < y.Length; i++)
                {
                    for (int j = 0; j < 3; j++)
                    {
                        F[i, j] = F[i, j + 1];
                    }
                    F[i, 3] = f.Get(i, x + h, y);
                }

                for (int j = 0; j < 1; j++)
                {
                    for (int i = 0; i < y.Length; i++)
                    {
                        y[i] = copy[i] + (9 * F[i, 3] + 19 * F[i, 2] - 5 * F[i, 1] + F[i, 0]) * h / 24;
                    }

                    for (int i = 0; i < y.Length; i++)
                    {
                        F[i, 3] = f.Get(i, x + h, y);
                    }
                }
            }
        }
    }
}
